#This file does two things:
#1) rich substitution in BLP data, ie. 
# elasticities, own and cross, MD,  and excess variances using BLP data
# this file will generate a nice avplot, and also a number of rich subs
# moments to be matched in our simulations.
#2) a CF on BLP data
#
rm(list=ls())
library(data.table)
library(HeadR)
library(doParallel)
library(doRNG)
library(SQUAREM)
library(fixest)
library(tictoc)
library(xtable)
library(latex2exp)
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_functions.R")
source("Rfiles/CF_MMX_functions.R")
source("Rfiles/CF_MMX_MLOG_functions.R")
library(Rcpp)
Rcpp::sourceCpp("Cppfiles/crossDelta.cpp")
#parallel processing set up for doParallel, doRNG
cl <- detectCores()-1
registerDoParallel(cores=cl)
getDoParWorkers()
MAXIT <- 500
seedn <- 777 # 666
sequences <- "normal" # alternatives: normal, halton or "sobol"
#
n.hh <- 100000 # c
U0 <- 1  #
n.reps <- 1 # 
#
###
# Targets 
###
source("Rfiles/CF_BLPdata_getdata.R") # get original data
BLP <- BLP[year==1990]  # from now on, we'll work with 1990 data only
BLP[,summary(100*s)]
#multiproduct firms
BLP[,m := 1:.N]
BLP[,.(num_models=uniqueN(m),num_firms=uniqueN(f))]
NM <- BLP[,.(num_models=uniqueN(m)),by=f]
NM[order(-num_models)][1:3][,sum(num_models)]/NM[,sum(num_models)] 
inside.share <- BLP[,sum(s)]
dom.share <- BLP[i==1,sum(s)]/inside.share
DT <- BLP[,.(y = sum(s)),by=f][order(-y)]
CR5.firms <- sum(DT$y[1:5])/inside.share
DT <- BLP[order(-s)]
CR5.models <- sum(DT$s[1:5])/inside.share
targets <- data.table(out.share = 1-inside.share,dom.share=dom.share,CR5.firms=CR5.firms,CR5.models=CR5.models)
sink(file="Tables/CF_BLPdata/BLPdata_targets.txt")
targets
sink()
#BLP stats table
BLPstats.tab <- data.table(round(100*targets,0),num_firms=uniqueN(BLP$f),num_models=uniqueN(BLP$m))
names(BLPstats.tab) <- c("Outside good share (\\%)","Domestic share (\\%)","Concentration (CR5 in \\%) firms","Concentration (CR5 in \\%) models","Number of firms", "Number of models")
BLPstats.tab <- data.table(BLPstats.var=names(BLPstats.tab),BLPstats.value=as.numeric(BLPstats.tab[1]))
BLPstats.tab[, order:=1:length(BLPstats.tab$BLPstats.var)]
names(b.means) <- c("Constant","HP/WT","Air con.","Miles/USD","Size")
b.tab <- data.table(var=names(b.means),b=b.means,sd=b.sds)
mu <- y.mn[year ==1990]$V2
SD.alpha_h = alpha*sqrt((exp(y.sd^2)-1)*exp(-2*mu+y.sd^2))   
b.tab <- rbind(b.tab,data.frame(var="Price",b=alpha,sd=round(SD.alpha_h,3)))
b.tab[, order:=1:length(b.tab$var)]
BLPstats.tab <- merge(b.tab,BLPstats.tab,by="order",all.x = TRUE) 
BLPstats.tab[,output:=texout(list(var,b,sd,BLPstats.var,BLPstats.value),3)]
sink(file="Tables/CF_BLPdata/BLPstats_tab.tex")
writeLines(BLPstats.tab$output)
sink()
#
###
# What is the average own elasticity? ====
###
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_getdata.R") # 
param99 <- fread("Data/BLP99/published_param_1999.csv") # BLP1999 parameters
param99[, xvar := st_right(param_name,"_")]
alpha99 <- param99[param_name=="alpha_price"]$pub_param
b.means99 <- param99[param_name %like% "demand"]$pub_param
names(b.means99) <- param99[param_name %like% "demand"]$xvar
b.sds99 <- param99[param_name %like% "sigma"]$pub_param
names(b.sds99) <- param99[param_name %like% "sigma"]$xvar
#
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] # 
BLP.meanownelas
#
###
# Compute own elas for a given variety ====
###
#
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_getdata.R") # get original data
#
BLP[,rank := frank(-s),by=year]
View(BLP[year==1990,.(m,p,s,rank)][order(-s)])
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] 
BLP.priceadjust <- BLP.output.m[, mean(p_m * (1-s_m))] 
#
# graph of elas vs price for many setting versions for chosen model (currently Volvo) 
shock = 0.25
model.chosen <- "VV24089" 
#
for(vsim in 1:3) {
  set.seed(seedn)
  data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=vsim) 
  BLP.output <- xi.BLP.mlog(dat=data.exog.mm,verbose=FALSE) 
  model.eq <- elas_price(p.change=0,dat=data.exog.mm,xi.dat=BLP.output,m.chosen=model.chosen)
  DT <- lapply(seq(from=-shock,to=shock,by=.01),elas_price,dat=data.exog.mm,xi.dat=BLP.output,m.chosen=model.chosen)
  assign(paste0("DT",vsim),rbindlist(DT))
  }
DT3[,epsoly := -BLP.meanownelas*(1-s) ]
saveRDS(DT1,"Data/RDS/DT1.rds")
saveRDS(DT2,"Data/RDS/DT2.rds")
saveRDS(DT3,"Data/RDS/DT3.rds")
DT1 <- readRDS("Data/RDS/DT1.rds")
DT2 <- readRDS("Data/RDS/DT2.rds")
DT3 <- readRDS("Data/RDS/DT3.rds")
#should saveRDS these DT1, DT2, DT3, because they are small but time-consuming to create
# Fig 2(a) (Fig 2(b) generated in CF_BLPdata_ysdloop.R)
brilliantblue =rgb(0,117,197,maxColorValue = 255)
my.darkgray <- adjustcolor("black",0.7)
#margins to use everywhere
pltmgn = c(4,4,1,1)
pdf(file=paste0("Graphs/CF_BLPdata/elas_price_settings_",model.chosen,".pdf"),height=4,width=4.5)
par(mar=pltmgn)  #bottom, left, top and right
yrng <- range(DT1$eps)
symbvals = c(-0.25,-0.13,0,0.13,0.25)
pvals <- DT1[p.change %in% symbvals]$p
#Mixed logit
DT3[,plot(p,eps,log="xy",col=my.darkgray,type="l",ylim=yrng,
          xlab="price (log scale)",ylab=TeX("Own price elasticity ($\\epsilon$, log scale)"))]
#points(pvals,DT3$eps[which(DT3$p %in% pvals)],pch=16,col=my.darkgray)
points(pvals,DT3[p %in% pvals]$eps,pch=16,col=my.darkgray)
#CES-MC
abline(h=-BLP.meanownelas,col="black")
#CES-OLY
DT3[,lines(p,epsoly,col="black",lty="dashed",lwd=2)]
#Beta-het
DT2[,lines(p,eps,col=my.darkgray)]
points(pvals,DT2[p %in% pvals]$eps,pch=15,col=my.darkgray)
#Logit
DT1[,lines(p,eps,col=my.darkgray)]
points(pvals,DT1[p %in% pvals]$eps,pch=2,col=my.darkgray,cex=1.75)
#legend for lines
legend("topleft",
       col= c(my.darkgray,my.darkgray,my.darkgray,"black","black"),
       pch= c(NA,NA,NA,NA,NA),
       lty=c("solid","solid","solid","solid","dashed"),
       legend=c("Mixed Logit","Logit",TeX("Logit with $\\beta$ het."),"CES-MC","CES-OLY"),
       lwd=1.5,cex=0.7
       )
#legend for symbols
legend("topleft",
       col= c(my.darkgray,my.darkgray,my.darkgray,NA,NA),
       pch= c(16,2,15,NA,NA),
       lty=c("solid","solid","solid","solid","dashed"),
       legend=c("Mixed Logit","Logit",TeX("Logit with $\\beta$ het."),"CES-MC","CES-OLY"),
       lwd=1.5,cex=0.7,
       pt.cex =c(1,1.5,1,NA,NA), bty = "n"
)
dev.off()

#
###
# Implement CF on BLP ====
###
#
source("Rfiles/CF_BLPdata_getdata.R") # get original data
year.gph <- 1990
tar <- 0 # original tariff (to be raised)
tar.increase <- 0.1
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix 
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#do a settings file only for the speed 
settings.var <- data.frame(speed =c(0.3,0.3,0.3))
#
# Run the calibration to get xi and mc
year.monte <- year.gph
objective <- "CF"
market <- 1
n.x <- 4
#
#Run sims for the 3 versions
#EHA:
set.seed(seedn)
CF_BLPdata.simulv1 <- CF_BLPdata.simul(vsim=1,cf.meth="EHA")
set.seed(seedn)
CF_BLPdata.simulv2 <- CF_BLPdata.simul(vsim=2,cf.meth="EHA")
set.seed(seedn)
CF_BLPdata.simulv3 <- CF_BLPdata.simul(vsim=3,cf.meth="EHA")
#
CF_BLPdata.simulallv <- rbind(CF_BLPdata.simulv1,CF_BLPdata.simulv2,CF_BLPdata.simulv3) 
saveRDS(CF_BLPdata.simulallv,"Tables/CF_BLPdata/CF_BLPdata_simul.rds")
#
#
#AHA:
set.seed(seedn)
CF_BLPdata.simulv1 <- CF_BLPdata.simul(vsim=1,cf.meth="AHA")
set.seed(seedn)
CF_BLPdata.simulv2 <- CF_BLPdata.simul(vsim=2,cf.meth="AHA")
set.seed(seedn)
CF_BLPdata.simulv3 <- CF_BLPdata.simul(vsim=3,cf.meth="AHA")
#
CF_BLPdata.simulallv <- rbind(CF_BLPdata.simulv1,CF_BLPdata.simulv2,CF_BLPdata.simulv3) 
setnames(CF_BLPdata.simulallv,old=c("chg_ces_qty","bias_qty"),new=c("chg_ces_qty_AHA","bias_qty_AHA"))
CF_BLPdata.simulallv <- CF_BLPdata.simulallv[,list(v,chg_ces_qty_AHA,bias_qty_AHA)]
saveRDS(CF_BLPdata.simulallv,"Tables/CF_BLPdata/CF_BLPdata_simul_AHA.rds")
#
###
# Implement Rich subs ====
###
#
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_getdata.R") # get original data
objective <- "richsubs"
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2) # this is the median income , not the average
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#
objective <- "CF" # if marginal costs are needed 
source("Rfiles/CF_BLPdata_xi.R") # extract the implied xi
#
source("Rfiles/CF_BLPdata_richsubs.R")  # cross-elasticities etc.


# Table for slides and paper####
#
CF_BLPdata.simulallv <- readRDS("Tables/CF_BLPdata/CF_BLPdata_simul.rds")
CF_BLPdata.simulallv_AHA <- readRDS("Tables/CF_BLPdata/CF_BLPdata_simul_AHA.rds")
CF_BLPdata.simulallv <- merge(CF_BLPdata.simulallv,CF_BLPdata.simulallv_AHA,by="v")
monte.BLPmoments <- readRDS("Tables/CF_BLPdata/BLPrichsubs_moments.rds")
CF_BLPdata.table <- merge(CF_BLPdata.simulallv,monte.BLPmoments,by="v")
CF_BLPdata.table[, v_name:="Mixed Logit"]
CF_BLPdata.table[v==1, v_name:="Logit"]
CF_BLPdata.table[v==2, v_name:=" $\\beta $ het."]
CF_BLPdata.table[v_name=="Mixed Logit", v_table:=1]
CF_BLPdata.table[v_name=="Logit", v_table:=2]
CF_BLPdata.table[v_name==" $\\beta $ het.",  v_table:=3]
#
CF_BLPdata.table <- CF_BLPdata.table[,list(v_table,v_name,meanownelas=-meanownelas,meansuperelas,pt,pte,pte1,
                                           chg_true_qty,chg_ces_qty,chg_ces_qty_AHA,bias_qty,bias_qty_AHA
                                           ,crosselasmaha,exvar,covprob)]
CF_BLPdata.table.DGP <- CF_BLPdata.table[,list(v_table,v_name,meanownelas,meansuperelas,crosselasmaha,covprob,exvar)]
CF_BLPdata.table.CFR <- CF_BLPdata.table[,list(v_table,v_name,chg_true_qty,chg_ces_qty,chg_ces_qty_AHA,pt,pte,pte1)]
#
#Outputs
setorder(CF_BLPdata.table, v_table)
CF_BLPdata.table[,c("v_table"):=NULL]
#
sink(file="Tables/CF_BLPdata/CF_BLPdata_simul_richsubs.tex")
print(xtable(CF_BLPdata.table,digits=c(0,rep(2,13),3)),sanitize.text.function = function(x){x}
      ,include.colnames = FALSE, include.rownames = FALSE, only.contents=TRUE, hline.after=NULL,comment=FALSE)
sink()
#
setorder(CF_BLPdata.table.DGP, v_table)
CF_BLPdata.table.DGP[,c("v_table"):=NULL]
#
sink(file="Tables/CF_BLPdata/CF_BLPdata_DGP.tex")
print(xtable(CF_BLPdata.table.DGP,digits=c(0,rep(2,5),3)),sanitize.text.function = function(x){x}
      ,include.colnames = FALSE, include.rownames = FALSE, only.contents=TRUE, hline.after=NULL,comment=FALSE)
sink()
#
setorder(CF_BLPdata.table.CFR, v_table)
CF_BLPdata.table.CFR[,c("v_table"):=NULL]
# a version with AHA
sink(file="Tables/CF_BLPdata/CF_BLPdata_CFR_AHA.tex")
print(xtable(CF_BLPdata.table.CFR,digits=c(0,rep(2,7))),sanitize.text.function = function(x){x}
      ,include.colnames = FALSE, include.rownames = FALSE, only.contents=TRUE, hline.after=NULL,comment=FALSE)
sink()
# a version without AHA
CF_BLPdata.table.CFR[,c("chg_ces_qty_AHA"):=NULL]
sink(file="Tables/CF_BLPdata/CF_BLPdata_CFR.tex")
print(xtable(CF_BLPdata.table.CFR,digits=c(0,rep(2,6))),sanitize.text.function = function(x){x}
      ,include.colnames = FALSE, include.rownames = FALSE, only.contents=TRUE, hline.after=NULL,comment=FALSE)
sink()





